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The theory underlying the Car-ParrineUo extended-lagrangian approach to ab initio molecular 
dynamics (CPMD) is reviewed and reexamined using 'heavy' ice as a test system. It is emphasized 
that the adiabatic decoupling in CPMD is not a decoupling of electronic orbitals from the ions but 
only a decoupling of a subset of the orbital vibrational modes from the rest of the necessarily-coupled 
system of orbitals and ions. Recent work ( J. Chem. Phys. 116 , 14 (2002) ) has pointed out that, 
due to the orbital-ion coupling that remains once adiabatic-decoupling has been achieved, a large 
value of the fictitious mass ^ can lead to systematic errors in the computed forces in CPMD. These 
errors are further investigated in the present work with a focus on those parts of these errors that are 
not corrected simply by rescaling the masses of the ions. It is suggested that any comparison of the 
efficiencies of Born-Oppenheimer molecular dynamics (BOMD) and CPMD should be performed 
at a similar level of accuracy. If accuracy is judged according to the average magnitude of the 
systematic errors in the computed forces, the efficiency of BOMD compares more favorably to that 
of CPMD than previous comparisons have suggested. 



I. INTRODUCTION 

Since its introduction by Car and Parrinello in 1985[H, 
the field of ab initio molecular dynamics (AIMD) has 
grown very rapidly. Although many of its initial suc- 
cesses involved its application to questions in condensed 
matter physics and materials science, it has now been ap- 
plied with success to a wide range of problems in diverse 
areas of science such as chemistry, biology and geophysics 
and contributed a great deal to our understanding of 
these fields. The central idea of the method is that the 
Kohn-Sham orbitals 2J that describe the electronic state 
of a system within density functional theoryQ (DFT) 
are evolved simultaneously with the ions as classical de- 
grees of freedom. In the original approach of Car and 
Parrinello, an inertial parameter is introduced in order 
to associate a Newtonian dynamics with the electronic 
orbitals. This parameter is an unphysical quantity that 
is commonly referred to as the fictitious mass, fi. While it 
is known that the introduction of this extra inertia into a 
system affects its dynamics, the extent to which dynam- 
ics are altered and the way in which they are altered are 
only beginning to be studied in detail 4, 5, 6]. 

In this paper some of the theoretical ideas that un- 
derlie the original, and still widely-used, Car-Parrinello 
approach to ab initio molecular dynamics are examined 
and the possible importance of the fictitious mass effects 
are investigated by studying the example of ice. In what 



follows, this approach to AIMD will be referred to as Car- 
Parrinello molecular dynamics (CPMD). In CPMD the 
equations of motion of the ions (with coordinates Rj and 
masses Mj) and the electronic orbitals \tpi) are derived 
from the classical lagrangian 



L = /i^(V',|^.) + i ^ MiRj - E[{^,}, {Rj}] 



(1) 
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and the motion is further constrained to the hyper- 
surface defined by the orbital orthogonality condition 
{tpi\ipj) = Siji^i, j^. Inspired by the ideas of Car 
and Parrinello, an alternative approach to AIMD was 
proposed^ in which the electronic orbitals are kept as 
close as possible to the instantaneous ground state using 
a combination of extrapolation from earlier times and 
explicit minimization. This general approach will be re- 
ferred to here as Born-Oppenhcimer molecular dynamics 
(BOMD). 

It has been known since its introduction that CPMD 
is an approximation to a fully converged BOMD simula- 
tion in which the electronic system is exactly in its ground 
state at every instant of the ionic motion, however CPMD 
is frequently preferred over BOMD because it is deemed 
computationally more efficient. In CPMD simulations, 
because the orbitals are treated as classical degrees of 
freedom on the same footing as the ions, there exists an 
energy (which includes the kinetic energy associated with 
the fictitious classical motion of the orbitals) that is per- 
fectly conserved as long as the time step that is used 
is small enough to adequately integrate the equations of 
motion. On the other hand, while the conserved quantity 
in BOMD is the physically meaningful total energy of the 
system of electrons and ions, its conservation in practice 
is always imperfect. There always exists a drift in the 
calculated total energy of a BOMD simulation because 
the electronic orbitals are never perfectly converged to 
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the ground state. While the magnitude of this drift can 
systematically be reduced by more fully converging the 
minimization of the electronic system at each time step, 
any such improvement in accuracy must always be bal- 
anced by the computational expense involved. The argu- 
ment has frequently been madeP.ITolllll|. that for many 
systems CPMD is computationally superior to BOMD 
because the computational overhead required to make 
BOMD conservative within acceptable limits is greater 
than that required to perform a perfectly conservative 
CPMD simulations with a value of /i chosen in line with 
conventional wisdom [sl 1^ fl^. Since BOMD can be as 
fast as, or faster than, CPMD depending on the size of 
the energy drift that is tolerated, an implicit assumption 
of this argument is that BOMD can only match the ac- 
curacy of CPMD when this drift of energy is very small. 
In this paper, this assumption is questioned by pointing 
out that, for commonly used values of the fictitious mass 
/i, there can exist large errors in the forces on the ions in 
a perfectly conservative and well-controlled CPMD sim- 
ulation. 

The emphasis of this paper is on the microscopic details 
of CPMD simulations, i.e. on the ability of the CPMD 
approach to evolving the electronic orbitals in time to 
produce (on average) the correct ground state forces on 
the ions. The important question of how deviations of 
forces from the correct forces relate to possible deviations 
in thermodynamic quantities from those that would be 
obtained using the correct forces is not addressed in the 
present work. 

Car and ParrincUo introduced their method in the 
context of coupling Kohn-Sham density functional 
theory[2,l3 to molecular dynamics, however, it has since 
been applied successfully to other electronic structure 
methods u^. It has also been applied in a purely classical 
context to evolve induced dipoles 19, 20] and fluctuating 
charges|22 in molecular dynamics simulations. None of 
these applications will be discussed in the present work 
but, clearly, many of the concepts discussed here have 
analogies within these methods. 



II. BACKGROUND 

It has been argued, and demonstrated for the case of 
silicon]^ ^|, that although there exist small instanta- 
neous deviations of the forces in a CPMD simulation from 
the ground state forces at the same ionic positions, these 
deviations average to close to zero on a femtosecond time 
scale and therefore do not result in serious errors in ther- 
modynamic properties^ Hit ll^ l2^ . The explanation 
proposed for this behaviour is that the electronic or- 
bitals perform small-amplitude ultra high frequency os- 
cillations around an equilibrium that slowly evolves as 
the ions move. If this equilibrium corresponds to the 
electronic ground state, and the amplitudes of the oscil- 
lations are small, the forces on the ions oscillate about the 
true ground state forces and errors in the forces average 



to zero on time scales that are relevant to ion dynamics. 

The fact that a dynamics is associated with electronic 
orbitals means that, according to the equipartition the- 
orem, the subsystems of orbitals and ions should equili- 
brate with one another and reach a common temperature. 
It has been observed, however, that for systems in which 
there is a substantial gap in the Kohn-Sham energy spec- 
trum between occupied and unoccupied electronic states, 
there is no systematic net transfer of kinetic energy be- 
tween orbitals and ions and therefore the orbitals remain 
at a temperature that is very low relative to the ionic 
temperature. This low temperature means that orbitals 
do not have the kinetic energy required to leave the en- 
ergy basin of the electronic ground state. 

The non-equilibration of orbitals and ions has been 
explained|l2j as being due to a gap in the vibrational fre- 
quency spectrum of the coupled orbital-ion system that 
separates the ultra high frequency orbital modes from all 
the lower frequency modes of the system. The ultra high 
frequency oscillations of the orbitals are "adiabatically 
decoupled" from the rest of the system. 

While the arguments presented above appear to work 
well for silicon^, j little work has been done to ver- 
ify that CPMD forces average to the ground state forces 
for other systems. Furthermore, some important details 
of the arguments of Pastore et al. have received little 
attention: Adiabatic decoupling is not a decoupling of 
orbitals from ions. Rather, it is only a decoupling of 
the quasi-independent orbital modes that consist of ul- 
tra high frequency oscillations of the orbitals from the 
lower frequency modes of the system. The lower fre- 
quency modes of the system generally involve both ions 
and orbitals. The slow (ionic time scale) component of 
the orbital motion is frequently cornpletely neglected in 
discussions of the CPMD method |5lllll|. however, and it 
is often stated that the lowest frequency motion in the 
orbital subsystem is approximately related to the Kohn- 
Sham energy gap of the system Eg through |^ ITTl IT^ 



(2) 



. In an insulator, this frequency is very high compared to 
typical ionic frequencies. Clearly, this view is not com- 
patible with the electronic orbitals remaining at or near 
the electronic ground state because the ground state by 
definition evolves on ionic time scales. 

The slow component of orbital motion due to the evo- 
lution of the ground state is always present and, as a re- 
sult, there is always a continuous exchange of energy and 
momentum between orbitals and ions. Furthermore, for 
many systems, this ionic time scale component of the or- 
bital motion has been found to be appreciable pi [fl Esjl . 
It will be demonstrated in section of the present work 
that this ionic-time scale contribution to the motion of 
the orbitals actually accounts for the vast majority of 
their classical "fictitious" kinetic energy in the simula- 
tions of ice reported here. 
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The argument of Pastore et al.^^ that forces in 
CPMD oscillate about the true ground state forces relies 
on the assumption that the ultra high frequency oscilla- 
tions of the orbitals are about the electronic ground state, 
however, for any finite /i it will be shown in section Fill Al 
that this is not the case. As a consequence of the ionic 
time scale motion of the electronic orbitals in CPMD, 
the orbitals do not oscillate around the ground state but 
about average values that are displaced from the ground 
state. Furthermore, although the resulting errors in the 
forces are small for materials in which valence electrons 
are delocalized and therefore have a low quantum kinetic 
energy (such as silicon), for other materials the forces 
can deviate strongly from their ground state values if 
commonly used values of the fictitious mass parameter 
are usedjj], however, these errors can systematically be 
reduced by reducing the magnitude of /x. 

It has previously been proposed that an effect of the 
orbital-ion cou plin g is simply to rescale the masses of 
the ionsj3, ITM Il6l|. an effect that does not alter ther- 
modynamic properties, however, theoretically, the error 
induced by orbital-ion coupling only reduces to a mass 
correction in the limit in which ions do not interact with 
one another, i.e. in the limit of an infinitely dilute gas|l5|. 
For condensed systems most, but not all, of the error can 
be corrected by simply rescaling the masses of the ions. 
One goal of the present paper is to review, to extend, and 
to help clarify some of the analysis of reference 0. An- 
other goal is to stimulate further investigation by using 
'heavy' ice as a test system to highlight the importance 
of the effects described in reference with particular em- 
phasis on the part of the errors in the forces that is nei- 
ther oscillatory on a very short time scale, or equivalent 
to a simple renormalization of the masses of the ions. It 
is shown that, for commonly used values of /i, the mag- 
nitude of this part of the force errors is much larger than 
errors in the forces in typical BOMD simulations. 

Water is arguably the most important system studied 
with ab initio molecular dynamics. CPMD has been ap- 
plied extensively to the study of liquid water [sl 1^. ITol Il7| 
and to the study of biological systems where water is 
present. The water molecules in ice have much the same 
properties as they do in liquid water, but ice has a larger 
electronic energy gap making it easier to simulate with 
CPMD. The goal here is to investigate the accuracy of the 
bare Car-Parrinello method without any further sources 
of error. In simulations of liquid water the smaller en- 
ergy gap between occupied and unoccupied states and 
the high frequency of the 0-H stretch vibrational mode 
can result in a continuous drift of kinetic energy into 
the orbital subsystem that would further complicate our 
analysis 0. These problems are avoided here by studying 
'heavy' ice rather than water. The larger energy gap of 
ice and the lower frequency of the 0-D stretch in D2O 
ensures that no such coupling occurs in the simulations 
reported here. 

The fact that the water molecules are confined to a lat- 
tice means that large deviations of the electronic struc- 



ture of the individual molecules from their average elec- 
tronic structure are less likely to occur in ice than in liq- 
uid water. As mentioned above, a rigid electronic struc- 
ture means that errors associated with fi only affect dy- 
namical properties whereas changes in the local electronic 
structure due to interactions between molecules may lead 
to errors in thermodynamic properties. For these rea- 
sons, the errors in CPMD reported here for ice are taken 
to be indicative of similar or more serious problems in 
liquid water. 

The fictitious mass effects described here and in refer- 
ence should not be confused with those described by 
Grossman et al. Q . The work of Grossman et al. pointed 
out some problems that occurred in simulations of water 
due to the adiabatic decoupling condition breaking down 
and energy passing continuously into the ultra high fre- 
quency orbital modes from the lower frequency modes. 
The present work is only concerned with fictitious mass 
dependent errors that are present in a well-performed 
Car-Parrinello simulation where the adiabatic decoupling 
condition is perfectly maintained and therefore the ultra 
high frequency orbital modes are energetically isolated. 



III. THEORY 

In this section the dynamics of the orbitals and the ions 
that result from the Car-Parrinello lagrangian of equation 
^are analysed in order to gain a better understanding of 
the different contributions to the errors in the forces on 
the ions. 

The equations of motion of the orbitals and ions that 
are derived from equation ^ are 



Miiif 



SE[{^P,},{-Ri}] 



dE[{^,},{Ki}] 



dRf 



(3) 



(4) 



where the orbital orthonormality condition is implicit 
in the functional derivatives 5/6ilj*{r) and all functional 
derivatives throughout this paper. The notation d/dRj 
indicates that a partial derivative is performed with re- 
spect to Rf with all other ion Cartesian coordinates fixed 
and with the orbitals ipi fixed. In the remainder of the 
paper it will also be necessary to use the notation V" 
to denote a partial derivative with respect to Rf with 
all other ion Cartesian coordinates fixed but not with the 
orbitals fixed. 

Equation ^ produces the correct Born-Oppenheimer 
dynamics if the set of orbitals {ipi} at which the deriva- 



tives 



dEl{i,,}.,{-R,}] 
dRf 



are evaluated is the set of ground state 
orbitals {?Af'*'}, i.e. the set of orbitals that minimize 
the Kohn-Sham energy functional £'[{?/'i}, {R-/}] for fixed 
ionic positions {R/}- The quality of the forces on the ions 
therefore depends on the ability of equation|21to produce 
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a dynamics for the orbitals that keeps them close to this 
instantaneous ground state as the ions move. In section 
liii Al therefore, the dynamics of the orbitals are analysed 
in the regime of small deviations from the ground state 
in order to gain a better understanding of the different 
contributions to the errors in the forces on the ions. In 
section IIII ti\ the effects of the orbital dynamics on the 
forces on the ions are discussed. 



A. The dynamics of the orbitals 

As stated previously, a common view of Car-Parrinello 
dynamics is that the Car-Parrinello orbitals {i/ji} oscil- 
late rapidly around the electronic ground state so that 
deviations from the ground state average to close to zero 
on the longer time scales that are relevant to ionic dy- 
namics. In this section it is demonstrated that this is not 
the case if /i is large, and an expression for the average 
value of the electronic orbitals is derived. 

For a given set of ionic coordinates {R/}, the elec- 
tronic ground state is well defined. It is therefore con- 
venient to introduce the quantity \6^pi) = IV'i) — l^-'f''^')? 
i.e. the deviation of the ith Car-Parrinello orbital IV^i) 
from its ground state value. The state of the system 
at a given instant of a Car-Parrinello dynamics can be 
completely specified by specifying the values of the dy- 
namical variables {R./}, {Sipi} and their time derivatives 
{III} , {6ipi} at that instant. In what follows, the dy- 
namics of the {Sipi} will be analyzed. An assumption 
that will be made throughout the analysis is that the 
{STpi} are small enough that certain ^-dependent quan- 
tities are well described by truncations at linear order 
in Sip of Taylor expansions about the electronic ground 
state. 

A Taylor expansion of 5E/5%}j*{y) about the ground 
state gives 



5E 



SE 



E 



5^E 



S^E 



g.s. 



dv' + 0{5^?) (5) 



where the first term on the right hand side vanishes by 
definition of the ground state. 

The second time derivative of may be written as 



(6) 



where 



I,a 



I,a,J,/3 



(7) 



Clearly, ipf''^' (r) is finite unless the ions are not moving. 
At a given instant in time V'f ^ (r) depends on the ionic 



positions {R/}, the velocities {R/}, and on the orbitals 
{S^pi}. The dependence of '(/'f '^ (r) on the {S^pi} comes 
in via the dependence of the ionic accelerations on the 
{Sipi} through equation^ However, equation^ may be 
rewritten as 



Rf 



+ 



_! d_ 

Mi dRJ 

SE 



E 



g.s. 




M 



-V^E 



OiSiP^) 



(8) 



from which it is clear that the dependence of R", and 
therefore of tpf^'ir), on {S'lpi} is at higher than linear 
order. To first order in {Sipi}, therefore, equations 13 [H 
and El can be combined to give 



S'^E 



srA-r') 



+ 



S'^E 



g.s. 



(9) 



By assuming that the {■0i} are either real, or are complex 
but with phase factors that do not vary as the system 
evolves, this equation can be evaluated as 



where 



j K{iv,jv')5i:,{v')dv' -i^r-{v) (10) 



5v{v) 



+ 2 



5n{v') 



g.s. 



(11) 



g.s. 



H is the Kohn-Sham single particle Hamilitonian, v{y) = 
va (r) -I- wxc (r) is the sum of the Hartree and exchange- 
correlation potentials, and fj is the occupation number of 
the j*^ orbital. It is clear that K{i r,j r') is independent 
of the {Sipi} and therefore evolves on ionic timescales. 

If the {Sipi} are sufficiently small, their motion is gov- 
erned by eauation II 01 which describes driven coupled os- 
cillations of the {Stpi}, where the driving term is given by 
the acceleration of the ground state {V'f "'}• By means 
of a suitable unitary transformation, K{i r,j r') can be 
diagonalized and eauation 1101 recast into a simpler form 
involving transformations of the functions {dipi{r)} and 
{■0f '* (r)}, however, for the purposes of the present dis- 
cussion it is sufficient to notice that, by reducing the 
value of the free parameter fj,, one can make the magni- 
tude of K/ ji arbitrarily large thereby increasing the fre- 
quencies of the oscillations to "ultra high" values, while 
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the driving term ipf''^'(r) remains constant. It is assumed 
here that a value of /i has been chosen such that there 
exist three distinct time scales rg r/ <C tl in the 
Car-Parrinello dynamics. T5 is a short time scale that 
is comparable to the period of the ultra high frequency 
oscillation of the orbitals, is the long time scale on 
which the ions move, and t/ is an intermediate time scale. 
When averaged over a time scale of order ^ r/, equation 
IIUI becomes 



For larger values of /i, or if the particulars of a simula- 
tion are such that ultra high frequency oscillations of the 
{STpi} have a large amplitude, contributions to the errors 
in the forces that are of higher than linear order may 
play a role, even if the separation of time scales rg tl 
can still be maintained, however, such effects will not be 
discussed in the present work. 

B. The forces on the ions 



SMr) = --Y. I K{ir,jr')5^^{v')dv' -i:f'-{v) (12) 

where the averaging only significantly affects (r) and 
Sipi{r) because both K and i/'f'* (r) vary on the much 
longer time scale tl- If the orbitals are to remain close 

to the ground state, then Stpi « 0, Vi, because the short 
time scale dynamics of the orbitals should be oscillatory 
rather than accelerating systematically in one direction. 
By inverting equation I12L the average deviation of the 
wavefunctions from their ground state values can there- 
fore be written as 



K """(j r', i r)'0f (r)(ir 



(13) 



S^j{r') is, in general, non-zero and dependent only on 
ionic positions and velocities. Therefore, the electronic 
orbitals do not average to their ground state values on 
time scales much longer than the time scale of their ul- 
tra high frequency oscillations but shorter than the time 
scale of ionic motion. On the other hand, S^j{r') de- 
pends linearly on the fictitious mass /i and therefore can 
be made arbitrarily small by reducing its value. Further- 
more, as pointed out by Pastore et a/.0], by decreasing 
the value of fj, the frequencies of the oscillations of the 
Sipi about their slowly evolving average values can be 
made so large that virtually no energy is transferred to 
them from the lower (ionic) frequency modes of the or- 
bitals and ions. If the amplitudes of these oscillations 
are initially small, therefore, they will remain small un- 
less the system is artificially perturbed (by, for example, 
discontinuously changing the velocities of the ions). 

To summarize this discussion of the orbital dynam- 
ics: it has just been shown that for small values of the 
fictitious mass, /i, the dynamics of the Car-Parrinello or- 
bitals can be described as consisting of ultra high fre- 
quency oscillations about average values that evolve on 
ionic timescales and that are given by equation 1131 For 
sufficiently small values of /i ( and {Sipi}) the deviations 
of the forces on the ions from their values at the elec- 
tronic ground state can therefore be separated into an 
ultra high frequency oscillatory part that averages to zero 
on the time scale relevant to ion dynamics and a system- 
atic part that does not. The systematic part of the errors 
in the forces on the ions is the main focus of the present 
work. 



In this section, the impact of the orbital dynamics dis- 
cussed in the previous section on the forces on the ions 
are discussed. Of primary concern is the systematic part 
of the error in the forces due to the displacement of the 
average values of the orbitals from the ground state, 6ip. 
This contribution is now derived. 
Equation 0] can be written as 
dE 



dRf 



E 



5E di!*(v) 6E d^Jiir) 



(5?/'*(r) dRf Si^^ir) dRf 



dr 



Substitution of equation |21 yields 



F? = -VfE-Y,f, J 



^,(r)V?^A*(r)+^A*(r)V?V^(r 



(14) 



dr 



(15) 

since Vj'!/'i(r) = dil)i{r) / dRf . The first term on the 
right-hand side of eauation llSl is now expanded in a Tay- 
lor series about the ground state : 



VfE = -Vf< E 



6E 



E 



g.s. 



SE 



^BOi 



g.s. 

0- 



dY + 0{5^1) 



g.s. 



(16) 



where F^q^ is the a-th Cartesian component of the Born- 
Oppcnheimer force on atom /. As pointed out in section 

nil Al on time scales relevant to ionic motion V'i(r) — 

S^i{r) + ipf''^' {r) — ipf''''{r), therefore, equations 171 IT^ 
and 1161 can be combined to give the error in the Car- 
Parrinello force on timescales ( > t/) relevant for ionic 
motion aspsj 



F\ 



BOi 



?Vr-^-(r)V>f--(r) 



+ E 

J,P,K,'i 



^-•(r)VlV^Vf^-(r) 



dr 



(17) 
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This expression is vahd in the hmit of small {^V'l} and 
depends only on ionic positions and velocities. It has no 
first order dependence on the {Sipi}. 

This systematic error in the Car-Parrinello forces arises 
from the fact that the electronic ground state evolves 
on ionic time scales. As ions move, they must push 
the orbitals around. The orbitals carry inertia and so 
the ions must impart momentum to them in order for 
them to follow the ground state. This means that the 
ions lose momentum to the orbitals, i.e. the ions en- 
counter a resistance to their motion. It is as though the 
ions move through a viscous, inhomogcncous and time- 
varying medium. 

The right hand side of equation 1171 depends on the 
details of the electronic ground state, as well as the ve- 
locities of the ions, and so it is difficult to make general 
statements about how AF affects the dynamics and ther- 
modynamics of a simulation. In the absence of further 
information there is no reason to expect either dynam- 
ical properties or thermodynamic properties to be cor- 
rect if fi is large enough to make the error of equation 
1171 significant. However, in the limit in which individual 
ions are spherically-symmetric and do not interact with 
one another, this error has been shown to reduce to the 
form AF^ — —AMRf, where AM/ is a constant that is 
proportional to the total quantum kinetic energy of the 
electronic states on ion In other words, if ions 

don't interact with one another (i.e. in the limit of the 
infinitely dilute gasfl^) the systematic error in the force 
due to the displacement of the average values of the or- 
bitals from the ground state, {Sipi}, has the same effect 
as simply increasing the masses of the ions by an amount 
AM/. This result had been noticed much earlier jl^, and 
described in reference^^ The more general result of ref- 
erence was also independently found by BlochlTE'l. An 
equation of motion for the ions that is an approxima- 
tion to the equation of motion of the ions in CPMD is 
therefore given bv|l.l25l| 

(M/ + AM,)Rf . = FSo. (18) 

The degree to which equation^] describes the motion of 
the ions in CPMD depends on the degree to which the 
electronic structure of the condensed system in question 
can be described in terms of a linear superposition of 
tightly bound ionic orbitals that are transported rigidly 
(i.e. without changing shape in response to their environ- 
ment) as the ions move. This approximation to CPMD 
is known as the "rigid-ion" approximation!^]. Equation 
^|can be understood as arising from the fact that such 
tightly bound orbitals have inertia (oc /Lt) and that this 
adds to the total inertia of the ions: the ions have to 
carry the "heavy" orbitals. 

In a dynamics described by equation I18L i.e. a Born- 
Oppenheimer dynamics but with increased ionic masses, 
thermodynamic properties are unchanged. This is only 
the case, however, if thermodynamic quantities such as 
temperature and thermal pressure are calculated using 



these increased values of the masses. In reference this 
was demonstrated for the case of MgO: The true temper- 
ature of the oxygen subsystem was that calculated using 
a mass for the oxygen ions of Mq + AMq, where Mq was 
the bare oxygen mass of 16 a.m.u. and AMq a correc- 
tion that accounted for the inertia of the orbitals moving 
rigidly with the oxygen ions. Car-Parrinello simulations 
that neglect the correction to the temperature due to in- 
creased ionic masses have effectively been performed at 
temperatures that are higher than those reported. Fur- 
thermore, this correction is important if more than one 
ionic thermostat is used in a simulation. For example, if 
mass corrections are not accounted for in the calculation 
of temperature, and if a different thermostat is applied 
to each ionic species, the result can be that each species 
is effectively maintained at a different temperature. The 
definition of temperature may be particularly important 
if so-called "massive thermostating" is employed 26], in 
which each degree of freedom is coupled to its own ther- 
mostat. 

In a system in which there is only one atomic species 
and every atom is in an identical chemical environment 
and if the rigid-ion approximation works well, dynamical 
properties can be corrected simply by rescaling time by a 
factor MjiM -I- AM). If the rigid-ion approximation 
works, but there is more than one atomic species, dy- 
namical properties can only be made correct by rescaling 
the mass of each atom a 'priori. In other words, the simu- 
lation is performed with a nominal mass for each ion 1 of 
Ml — AM/ so that the mass of the ion once 'dressed' by 
the heavy Car-Parrinello orbitals is the correct mass M/. 
If this is not done, large errors in dynamical properties 
can result [3. 

Strictly speaking, the rigid ion approximation is not 
applicable when ions interact with one another. In other 
words, when there exist electronic orbitals that depend 
on the positions of more than one ion, i.e. for which 
V/(/) ^ and Vj0 ^ for some ions 1 ^ J . Of course, 
this is generally the case for all condensed matter and 
therefore for all systems of interest. What this means 
is that the dynamics and thermodynamics of CPMD al- 
ways differ from BOMD. Fortunately, for all condensed 
systems tested to date, most of the error of equation 1171 
can be corrected by applying mass corrections to the ions. 
The concern that remains, however, is that even though 
it is much smaller, the part of the error that remains once 
mass corrections have been applied may not be negligible. 

A further difficulty associated with correcting the 
masses of the ions is that the quantum kinetic energy 
of the electronic states on a given atomic species is not a 
well defined quantity. This is because each orbital can- 
not, in general, be associated with only one atom. A 
method must therefore be devised to approximate the 
mass corrections. In the simulations of D2O presented 
here, the mass corrections that are used are those that 
minimize the errors in the forces on the ions. 

In order to differentiate between the part of AF that 
amounts to a simple mass correction, and the part that 
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FIG. 1: Some samples of force components on oxygen ions 
along a 4 f s segment of CPMD trajectory 1 ps after temper- 
ature has been adjusted using velocity rescaling. The Car- 
Parrinello force and the Car-Parrinello force once corrected 
with the best constant mass correction is compared to the 
ground state (Born-Oppenheimer) force, FbOi the same 
ionic positions. 



remains once this mass correction has been applied, the 
quantities AF™ and AF'' are introduced where AF™ 
—AMR and AF'' = AF — AF™ is the remaining error. 
The expression for AF in equation 1171 has been derived 
by taking an average over the ultra high frequency orbital 
oscillations, and therefore a more complete expression for 
the Car-Parrinello force is 

F/ = Fso.({R/}) + AF/({R,},{R,}) 
+ (5F,({<5VJ) + OW) 

« Fso. +AF7^ + AFJ + ^F/ (19) 

where (5F/ is the oscillatory part of the force which, for 
small enough /i, averages to zero on time scales r/ 3> r^. 

Figure n demonstrates that, indeed, the difference be- 
tween the Born-Oppenheimer and Car-Parrinello forces 
consists of both a systematic part AF, and an oscillatory 
part (5F. The systematic part appears to be mostly cor- 
rected by the application of mass corrections to the ions. 
The forces in figure Q are those from a Car-Parrinello 
simulation in which velocity rescaling was employed for 
0.5 ps to bring the temperature from T ~ 120-?^ to 
T ~ 220 K. The velocities of the ions were adjusted 
only 4 times, in total, during this 0.5 ps . The system 
was then equilibrated for one picosecond after which the 
forces were examined along the 4 fs segment of trajec- 
tory and compared to the forces at the same ionic po- 
sitions, but with the orbitals in their ground state, i.e. 
the Born-Oppenheimer forces F bOi ■ The amplitude of 
the oscillations, SF, in figure ^ appears worryingly large. 
It will be shown in sections llVI and IVll however, that if 



care is taken to ensure that the velocities of the ions are 
not changed discontinuously, the large amplitude of these 
oscillations disappears and the contribution of SF to the 
instantaneous error in the Car-Parrinello force can be 
neglected. At higher temperatures, however, large am- 
plitude oscillations may always be present, as was the 
case in simulations of MgO in reference '3> 

In section lvTI the systematic part of the error AF, will 
be examined. It has been assumed in the past 10] that 
AF w AF'" and that AF'' can be neglected. Therefore, 
in section IVTl the validity of this assumption will be in- 
vestigated. The degree to which AF can be eliminated 
by applying mass corrections to the ions will be tested. 
In other words, the magnitude of AF'', which has the po- 
tential to alter thermodynamic properties and for which 
no correction yet exists, will be examined. 



IV. SIMULATION DETAILS 

In the CPMD simulations of D2O reported here norm- 
conserving pseudopotentials have been used to describe 
oxveenp^ and deuterium|23 and the valence electronic 
orbitals were expanded in plane waves with a maximum 
energy of 70 Rydberg. Simulations were performed on 
a 24 X 24 X 12 a.u. simulation box containing 32 D2O 
molecules. Only the F— point was used to sample the 
Brillouin zone. A gradient-corrected functional was used 
to treat the effects of exchange and correlation 30J . In a 
preliminary CPMD simulation of liquid water mass cor- 
rections of AMo = 6.77/i a.u. and AMu = 0.213/j, a.u. 
, were found for the oxygen and deuterium ions, respec- 
tively. These mass corrections were found by calculating 
the ground state forces at selected points along a trajec- 
tory and by minimizing the error in the CPMD forces 
relative to these ground state forces. These mass cor- 
rections were then used in all of the CPMD simulations 
of ice reported here by reducing the values of the ionic 
masses used in calculating the acceleration from New- 
ton's equation of motion i.e. 



Rf 



1 



dE[{ip,},{-Ri}] 



Mi - AM/ 



dRf 



(20) 



The values of the fictitious mass used in the simulations 
reported here were 900 a.u. and, for comparison, 100 
a.u. The equations of motion for electrons and ions were 
integrated using a molecular dynamics time step of 6 a.u. 
for the /i — 900 a.u. simulations and 2 a.u. for the /i = 
100 a.u. simulations. No mass preconditioning scheme 
was used in any of the simulations reported here. 

Much care has been taken to minimize errors in the 
simulations reported here. For example, as pointed out 
by Remler and Madden [sif. and discussed in section lTlI Bl 
it is important that the velocities of the ions are not 
changed discontinuously as this can give a sudden kick 
to the electronic orbitals that results in large-amplitude 
ultra high frequency oscillations of the orbitals around 
their equilibria. Great care has therefore been taken to 
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FIG. 2: A schematic illustrating the sequence of Car-Parrinello molecular dynamics simulations that were performed. 



eliminate this source of error from the simulations re- 
ported here with the result (evident in figures [7| and |S1 of 
section lyijl that instantaneous errors in the forces due to 
ultra high frequency oscillations of the orbitals (i.e. SF 
) are small enough relative to the systematic errors of 
equation El that they may safely be neglected. 

The input coordinates for the CPMD were obtained by 
performing a very long molecular dynamics simulation of 
ice at low temperature (~ 100 K ) using an ab initio 
parametrized polarizable atomistic potential of the same 
form as that constructed for silica in reference Is^. This 
potential does not provide a very realistic description of 
water but it was deemed preferable to using randomized 
initial coordinates. From these initial conditions a se- 
quence of CPMD simulations were performed as shown 
schematically in figure |21 CPMD simulations were be- 
gun with the electrons in their ground state and the ions 
at zero velocity. A fictitious mass oi fi = 900 a.u was 
used and therefore the ionic masses of oxygen and deu- 
terium were rescaled to Mq — AMq — 12.659 a.m.u and 
Md — AMu = 1.895 a.m.u, respectively and temperature 
was calculated using the true ionic masses of Mq — 16 
a.m.u and M]j = 2 a.m.u. After half a picosecond of sim- 
ulation in the microcanonical ensemble, the ions were at 
a temperature of 120 K. At this point two separate con- 
tinuations of the simulation were performed. The first 
continuation was performed in order to demonstrate the 
different contributions to the error in the forces on the 
ions and the effect of velocity rescaling on these forces, 
and has been discussed in section IIII Bl I n the second 
continuation, a Nose-Hoover thermostat 3^ was attached 
to the ions and the temperature was smoothly increased 
to approximately 250 K during a further half picosecond 
of simulation. The system was then equilibrated without 
a thermostat for 3.5 picoseconds. As shown in figurelHlthe 
temperatures of the subsystems of oxygen and deuterium 
ions were calculated and compared with and without the 
use of the corrections to the masses of the ions in the 
definition of temperature. As can be seen, when mass 




1.5 2 
Time [ps] 

FIG. 3: The temperatures of the oxygen and deuterium 
subsystems during the 3.5 ps equilibration simulation with 
fj, = 900 a.u. In (a) the temperatures have been calculated 
using masses for the ions that have been corrected to account 
for the extra inertia due to the weight of the Car-Parrinello 
orbitals (in this case 16 a.m.u and 2 a.m.u. due to the a 
priori rescaling of the masses). The oxygen and deuterium 
subsystems appear to be out of thermal equilibrium. In (b) 
the temperatures have been calculated in the standard way 
using the bare ionic masses (12.659 a.m.u and 1.895 a.m.u). 
The oxygen and deuterium subsystems appear perfectly equi- 
librated. 



corrections are used the oxygen and deuterium subsys- 
tems are at different temperatures indicating that the 
system is not well equilibrated. If mass corrections are 
not used, the system appears perfectly equilibrated. At 
this point, the simulation was again continued in two sep- 
arate simulations. In one of these simulations the value 
of /I = 900 a.u. was used as before. In another simula- 
tion the fictitious mass was changed to fi — 100 a.u. The 
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FIG. 4: The temperatures of the oxygen and deuterium sub- 
systems calculated using the corrected ionic masses (A/-I-AA/) 
in the two continuations of the simulation shown in figure |3 
In (a) a fictitious mass of 100 a.u. has been used and thermal 
equilibration between the subsystems is observed. In (b) a 
fictitious mass of 900 a.u. is used and the subsystems remain 
out of thermal equilibrium. 



rescaled masses of the ions in the simulation with /i = 100 
a.u. were 15.629 a.m.u. and 1.988 a.ni.u. If the rigid-ion 
approximation was perfectly applicable the fact that the 
ionic masses have been rescaled a priori would mean that 
the two simulations with different values of fj, should have 
identical dynamic and thermodynamic behaviour. 

It was found (figure 0)) that after more than a further 
2.5 ps the simulation with fi — 900 a.u. still showed no 
sign of thermalization according to the mass-corrected 
definition of temperature. However, the simulation with 
/i — 100 a.u. equilibrated quickly and the subsystems 
were at the same temperature after 1.5 — 2 ps (figure 0J. 
In order to compare the phonon spectra of ice using these 
two different values of for reasonably well equilibrated 
simulations at the same temperature it was decided to 
continue both simulations from the end of this initial 2.4 
ps run with /x = 100 a.u. For both — 100 a.u. and 
H = 900 a.u., a further 5.5 ps of simulation were carried 
out during which time the oxygen and deuterium subsys- 
tems remained at the same (mass-corrected) temperature 
in both simulations. For both values of /z, the velocity 
autocorrelation function was calculated on a time domain 
of 1.45 ps by averaging over this final 5.5 ps of simula- 
tion. These velocity autocorrelation functions were then 
fourier transformed to obtain the phonon power spectra. 

The reason for the inability of the ji = 900 a.u. sim- 
ulations to thermalize during the first 7 ps of simulation 
remains unclear. By using the dressed ionic masses in 
the calculation of temperature we are implicitly assuming 
that orbitals move rigidly with the ions and are unper- 
turbed by their environment. In a condensed system this 



is an approximation and deviations from the rigid-ion 
limit always occur. In the opposite limit to the rigid-ion 
approximation, i.e. in the limit of very weak orbital-ion 
coupling, the temperature of the ions should be calcu- 
lated with the bare ionic masses. For systems that are 
not perfectly ionic, therefore, the correct definition of 
temperature is unclear. While the results of section IVll 
suggest that a constant mass correction for the ions is not 
perfectly appropriate for D2O, the results of section El 
suggest that the rigid-ion approximation does a remark- 
ably good job of estimating the fictitious kinetic energy 
of the orbitals. In addition, the fact that during the final 
5.5 ps of simulation the oxygen and deuterium subsys- 
tems remained at the same mass-corrected temperature 
suggests that it is appropriate to calculated temperature 
using rescaled masses. 

The fact that reducing the value of /i appeared to facil- 
itate thermalization may indicate that the thermalization 
problem was an artifact of the fictitious mass fi, however, 
further work is required to clarify this issue. 

V. THE TIME SCALE OF ORBITAL MOTION 

A common view of CPMD is that, if a large enough 
gap exists between the energies of the occupied and un- 
occupied electronic states, the electronic orbitals move 
on time scales that are much faster than typical ionic 
time scales. For example, it has been suggested that the 
motion of the orbitals in CPMD can be approximately 
described as a superposition of oscillations whose fre- 
quencies are given by = (-^-^^7-^)^^^, where and ej 
are the Kohn-Sham energy eigenvalues of occupied and 
unoccupied electronic states, respectively fM Wi fllj. In a 
system with a substantial energy gap between occupied 
and unoccupied states, all these frequencies are very high 
compared to typical ionic frequencies. In this view of 
CPMD the lowest vibrational frequency present in the 
orbital subsystem is determined by the size of the energy 
gap and the fictitious mass. However, it should be clear 
from the work of Pastore et al. and sections lllland lllll of 
the present work that, if the orbitals are to remain close 
to the ground state and the ions are moving, their mo- 
tion must contain an ionic time scale component. Here 
it is shown that, in fact, it makes the dominant contribu- 
tion to the total orbital kinetic energy. The slow orbital 
motion results from the coupling between ions and the 
orbitals and it is precisely this coupling that leads to a 
systematic departure of the average forces in CPMD from 
the ground state forces. 

If the rigid ion approximation holds perfectly, then the 
part of the kinetic energy of the orbitals due to the evolu- 
tion of the electronic ground state can be obtained simply 
from the mass corrections and the velocities of the ions, 
i.e. 

E /^(^'l'^^) = ^ E ^MoRfR? + i E ^MoRfR? 

i ItO leD 
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FIG. 5: The fictitious kinetic energy (FKE) of the electronic 
orbitals as a function of time compared to the FKE predicted 
by assuming that orbitals are inert and rigidly follow the ions. 
The close agreement between the curves verifies that orbital 
motion occurs primarily on ionic time scales and that the 
contribution of the ultra-high frequency oscillations of the or- 
bitals about their instantaneous equilibria to the total FKE 
is tiny by comparison. 



(21) 

In figure |31the fictitious kinetic energy (FKE) as pre- 
dicted by the rigid ion model in equation I21l is compared 
to the FKE that is extracted from the CPMD simulation. 
The difference between them is also plotted and this is 
made up both of ultra high frequency oscillations of the 
orbitals about their instantaneous equilibria and the part 
of the ionic time scale evolution of the ground state that 
cannot be described within the rigid ion approximation. 
The very close agreement between the predicted FKE 
and the actual FKE demonstrates clearly that electronic 
orbitals move on ionic time scales and, furthermore, that 
this slow motion accounts for almost all of their ficti- 
tious kinetic energy. This also demonstrates that the 
lowest vibrational frequencies of the orbital subsystem 
are almost [33| independent of and are simply given 
by the lowest vibrational frequencies of the ionic subsys- 
tem. Some previous simulations have calculated the or- 
bital vibrational power spectrum by fourier transforming 
the orbital velocity autocorrelation function and have not 
detected such low frequencies [M IT^. However, in these 
simulations the orbital vibrations were analysed while the 
ions were kept stationary and therefore the contribution 
to the orbital motion from the evolution of the electronic 
ground state was not present. In reference IssI Herbert 
and Head-Gordon have analysed the orbital vibrations 
during ion dynamics, and figure 7 of that paper clearly 
shows a contribution to the orbital frequency spectrum 
that has a low (ionic) frequency and that appears almost 
independent of the fictitious mass. 



VI. FORCES 

The goal of first principles molecular dynamics is to 
make a direct connection between an accurate descrip- 



tion of the electrons and an accurate calculation of the 
physical property of interest. It is important to know 
that a calculated physical property that agrees well with 
experiment does so because the trajectory is realistic due 
to good forces calculated from a good treatment of the 
electronic structure. When there is a microscopic ba- 
sis for agreement with experiment the method can be 
applied with some confidence to situations or to prop- 
erties for which experimental data is unavailable. This 
point is stressed because agreement with experiment in 
molecular dynamics simulations can often occur due to 
a cancellation of errors. This has been observed both 
in first principles simulations 5] and in simulations us- 
ing empirical potentials where potentials that agree with 
experiment on many properties]^ have been found to 
provide a ver y p oor microscopic description of the inter- 
atomic forces |32j. 

It is obvious, therefore, that in order to proceed with- 
out recourse to empiricism or a posteriori experimental 
verification one should be able to depend on the accuracy 
of the computed forces. What is much less obvious is 
how one should judge the accuracy of forces. One way of 
quantifying the average magnitude of the departure of a 
set of forces {Fj} from some reference forces {F™^} is by 
computing the percentage root-mean-squared difference 
between the two sets of forces relative to the root-mean- 
squared force component, i.e. : 



A, 



= 100 X 



\Fr - 



cf||2 



El 



|^-cf|| 



(22) 



where the sum ^ is over a large number of ions / 
and over a large number of points v along the molecu- 
lar dynamics trajectory. In general, this is an extremely 
crude measure of the departure of the forces from the 
reference forces, however there is frequently no alterna- 
tive to using it. A number of classical potentials have 
been parametrized by minimizing this quantity while us- 
ing ground state density functional theory forces as the 
reference. It has been found that, for some simple oxides, 
values as low as Arms = 5 — 20% can be achieved and 
that, as long as Arms is sufficiently small (i.e. < 20% ), 
the accuracy of these force fields for many experimental 
quantities is well correlated with its value [33. l37| . Al- 
though it is crude and unsatisfactory, in the absence of an 
alternative quantitative general measure of the quality of 
forces in the presence of errors of unknown consequence, 
this quantity is used in the present work. 

In this section the forces on the ions in the CPMD 
simulations of ice are inspected along a segment of the 
trajectory. After approximately 1 ps of simulation the 
ground state forces were calculated along a 21 fs segment 
of the CPMD trajectory. These ground state forces were 
then used as the reference forces for calculating the root- 
mean-squared (r.m.s) relative error in the CPMD forces 
according to equation 1221 both before and after they had 
been corrected using the rigid-ion mass correction. In 
other words (using the notation introduced in section 
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nil Bl and assuming that SF" is negligible) the r.m.s. val- 
ues of AF" and AFJ " relative to the r.m.s. value of 
^BOi were computed. Before mass corrections were ap- 
plied it was found that, for oxygen ions, A,. amounted 
to 21.8%. In other words, the r.m.s. error in the forces 
on the oxygen ions was 21.8% of the r.m.s. oxygen ion 
force component. For the deuterons A^.m.s was 5.5%. 
The force on each atom is dominated by intra-molecular 
interactions which are relatively constant (in the refer- 
ence frame of the molecule) at room temperature and 
therefore are generally of much less interest to simulators 
than the inter-molecular interactions. For this reason, 
the most important quantity to examine is the net force 
on each water molecule. It was found that Ar.m.s for the 
net forces on the water molecules was 44.7%. The val- 
ues of the mass corrections for the oxygen and deuterium 
ions were varied so as to minimize Ar.m.s for the corrected 
forces. It was found that by applying mass corrections of 
AMo = 6.91/i = 3.41 a.m.u and AMd = 0.215fi = 0.11 
a.m.u the errors could be reduced substantially to 5.3% 
for oxygen, 1.4% for deuterium and 12.4% for the water 
molecules. These mass corrections are very close to those 
that were applied a priori based on preliminary tests of 
liquid water but in fact it was found that, for oxygen, 
a wide range of values of AMq gave very similar val- 
ues for the average error. The value of Ar.m.s. for the 
corrected oxygen and water molecule forces is plotted in 
figure El as a function of the mass correction for the oxy- 
gen ion. Ar.m.s. is found to be reasonably insensitive to 
the mass correction near the average optimal value be- 
cause the optimal mass correction varied from atom to 
atom or, for a given atom, it varied in time. By changing 
the value of AMq one improved the agreement of some of 
the forces with the ground state forces and disimproved 
others. It was also found that the mass correction that 
gave best agreement between ground state and corrected 
CPMD forces on a given atom was different, in general, 
for the different Cartesian components. In other words, 
the effective masses of the ions in this simulation were 
time-varying tensor quantities as is to be expected from 
equation El 

In figure [3 some sample force components on oxygen 
ions are shown to illustrate this fact. Figure[3(a) and fig- 
ure[3(b) are examples of oxygen ions which, over this seg- 
ment of trajectory, have different effective masses. Fig. 
13 (c) is an example of an oxygen ion force component for 
which the optimal mass correction varies considerably 
over this short length of trajectory. The large oscilla- 
tions of the Car-Parrinello forces of period ~ 1 fs that 
were visible in figure ^ are completely absent in figure 
demonstrating that the contribution, SF, of the ultra 
high frequency oscillations of the orbitals to the force is 
negligible due to the ions having been accelerated con- 
tinuously throughout this simulation. 

In figure|Slsome sample force components on molecules 
are shown and compared to the ground state force and 
the force once corrected using the average optimal mass 
corrections. 
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FIG. 6: The r.m.s relative error Ar.m.s in the mass-corrected 
forces on the oxygen ions (full line) and on the D2O molecules 
on the sample segment of the CPMD trajectory as a function 
of the magnitude of the mass correction, /i — 900 a.u. in this 
simulation. 
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FIG. 7: Some samples of force components on oxygen ions 
along a 21 fs segment of CPMD trajectory. The Car- 
Parrinello force (green) and the Car-Parrinello force once cor- 
rected with the best constant mass correction (blue) is com- 
pared to the ground state force at the same ionic positions. 
Also plotted in (a) and (b) are the Car-Parrinello forces once 
corrected using the optimal mass correction for that ion along 
this particular segment of trajectory. In (c) the optimal mass 
correction at three points along the trajectory are indicated. 
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FIG. 8: Some samples of force components on D2O molecules 
along a 21 fs segment of CPMD trajectory. The Car- 
Parrinello force (green) and the Car-Parrinello force once cor- 
rected with the best constant mass corrections (blue) is com- 
pared to the ground state force at the same ionic positions. 



It is clear from Fig. [7| and Fig. |S1 that the errors in 
CPMD forces due to the fictitious mass can not simply be 
described in terms of a constant correction to the masses 
of the ions. In other words, the {AFJ "} are quite large. 

To get some perspective on the level of errors in the 
forces, these errors are compared in magnitude to those 
at varying levels of convergence of the Kohn-Sham or- 
bitals to their ground state during a minimization us- 
ing the method of direct inversion in the iterative sub- 
space (DIIS) 38]. Fifteen equally-spaced atomic 'snap- 
shot' configurations were extracted from the final 5.5 pi- 
cosecond production run with /i ~ 900 a.u. and, start- 
ing with random wavefunctions, the Kohn-Sham energy 
was minimized. Convergence of this minimization was 
taken to mean that the difference in energy between suc- 
cessive electronic iterations dscf had fallen lower than 
10~^^ a.u. /atom. For each snapshot configuration all it- 
erations for which Sgcf > a.u. /atom were discarded. 
At all other iterations the percentage r.m.s errors in the 
forces Arms relative to the fully converged forces were 
calculated. A line was then fit to the dependence of 
logj^Q Arms on logj^Q Sscf ■ This was done separately for 
the forces on the oxygen ions, the deuterons and the wa- 
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FIG. 9: The r.m.s relative error Ar.m.s in the forces as a func- 
tion of the degree of convergence Ss.c.f of the total energy dur- 
ing minimizations to the ground state with DIIS for 15 'snap- 
shot' atomic configurations extracted from the CPMD simula- 
tions. Each circle corresponds to a single self-consistency step 
in the minimization of the wavefunctions for a single snapshot. 
The filled circles are the data from one full representative DIIS 
minimization and are included to illustrate the convergence 
characteristics for each snapshot. The dashed lines indicate 
the value of Ar.m.s. from the mass corrected and uncorrected 
forces in the CPMD simulations and the level of convergence 
of the total energy that this would correspond to. The grey 
line indicates the level of convergence that is typically used in 
BOMD simulations of water ( Ss.c.f = 10~^^ a.u. /atom ) and 
the value of Ar.m.s. that this corresponds to. 



ter molecules. The comparison with the errors in CPMD 
is quantitatively very similar in each case and therefore 
only the results for the oxygen ions are plotted in figure 
El Very similar results were also found when electronic 
minimization was performed using a conjugate-gradients 
technique and therefore it is assumed that, for the pur- 
pose of the present discussion, the dependence of Ar.m.s. 
on Ss.c.f could be assumed to be reasonably indepen- 
dent of the route taken to the electronic ground state. 
For comparison, the values of Arms for the forces from 
the CPMD simulations both before and after the appli- 
cation of the rigid-ion mass corrections are also shown 
in this figure El The results suggest that, for the mass 
corrected (uncorrected) forces, the degree of convergence 
required at each time step of a BOMD simulation in or- 
der to achieve an average error in the forces Ar.m.s. of 
the same magnitude is around Ss.c.f — 10~^ a.u. /atom 
{— 10~^ a.u. /atom ). This is many orders of magnitude 
larger than the level of convergence that is generally used 
for water ITol l39l l40l |. For example, in reference fioi 
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BOMD simulations were performed with a convergence 
of the energy at each step of Sg.c.f = 10~^^ a. u. /atom. It 
was reported that this simulation was about a factor of 
4 more expensive than a CPMD simulation that used a 
time step of 4 a.u. As illustrated in figure |51 the errors 
in the forces are about three orders of magnitude smaller 
than those in the CPMD simulation when Ss.c.f = 10^^^ 
a.u./atom^3|- A fair comparison between BOMD and 
CPMD should consider their relative efficiencies at the 
same level of accuracy or their relative accuracies at the 
same level of efficiency. In general, this comparison de- 
pends on the system under consideration and the level 
of accuracy (or efficiency) required. No attempt will be 
made here to study this question in detail. However, it is 
illuminating briefly to consider the relative efficiencies of 
BOMD and CPMD at the level of accuracy of a CPMD 
simulation with = 900 a.u., and at the level of accuracy 
of a BOMD simulations with Sg.c.f — 10~^^ a.u. /atom . 

It is assumed here that the efficiency of BOMD for ice 
is similar to that of water. Drawing on the results of Kuo 
et ai, it is therefore assumed that a BOMD simulation 
with Ss.c.f = 10~^^ is approximately three times slower 
than a CPMD simulation with /i = 900 a.u. and At ^ 6 
a.u. Assuming that errors in the forces scale linearly with 
/i, Ar.m.s. for the CPMD simulation can be reduced to 
the level of the highly converged BOMD simulation by 
reducing ii by about three orders of magnitudel4l'|. The 
time step required to integrate the equations of motion 
for the orbitals scales with the square-root of the ficti- 
tious mass, and therefore this reduction of requires a 
reduction of the time step (and therefore the speed of the 
simulation) by a factor of about 30. What this means is 
that, at this level of accuracy, BOMD is roughly an order 
of magnitude faster than CPMD. 

If a much lower accuracy is sufficient, such as the level 
of accuracy in a CPMD simulation with /i = 900 a.u., 
then a comparison with a BOMD simulation using a value 
of Ss.c.f that is five or six orders of magnitude larger 
is sufficient. At this level of convergence, the BOMD 
simulation of Kuo et al. would be faster, but because 
Sscf typically reduces by an order of magnitude with each 
one or two electronic iterations, the speed up may not 
be very great. It would depend on the average rate of 
convergence of the orbitals to the ground state and on the 
quality of the extrapolation of the orbitals from previous 
time steps. It seems likely that the speeds of CPMD and 
BOMD simulations would be much closer in this case. 
However, a BOMD simulation at this level of accuracy 
would suffer from discontinuous forces - a problem that 
is not present in CPMD. 

More detailed comparisons between CPMD and 
BOMD are clearly required. The point of the present 
discussion is simply that the accuracy (as judged by the 
quantity A„„s) of a reasonably well-converged BOMD 
simulation far exceeds that of a CPMD simulation using 
standard values (a few hundred atomic units) J2J of the 
fictitious mass. The argument that CPMD is more effi- 
cient than BOMdT^I is frequently based on comparisons 
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FIG. 10: The mass-corrected phonon power spectra for sim- 
ulations with ^ = 100 a.u. and ^ = 900 a.u. 

in which BOMD is held to a higher standard of accuracy 
than CPMD. 



VII. PHONON POWER SPECTRUM 

As described in section IIVI the phonon power spectra 
have been calculated for the simulations with /i = 100 
a.u. and /i = 900 a.u. The results are plotted in figure 
[TUl Different rescaled masses have been used in these two 
simulations and so the coincidence of the main peaks in 
these power spectra clearly vindicates the use of correc- 
tions to the masses of the ions. For MgO, it has already 
been demonstrated that the power spectrum is strongly 
altered by the fictitious mass if mass corrections are not 
applied 4]. In the present work only the corrected power 
spectra are presented, but it should be clear that, if the 
masses of the ions had not been rescaled a priori, the 
power spectrum of the simulation with /i — 900 a.u., in 
particular, would have been very different. Overall, the 
agreement between the two power spectra is reasonably 
good despite the fact that, as discussed in section IVII 
the masses of the oxygen ions vary considerably during 
the dynamics and that the mass correction that was ap- 
plied was simply a uniform constant correction to account 
for the average value of the extra inertia of the ions due 
to the fictitious mass. It may be that many time- or 
space-averaged properties are insensitive to /i as long as 
average constant mass corrections are applied to the ions. 
For example, the work of Grossman, Schwegler et aL[^Q 
suggests that once adiabatic decoupling is achieved, the 
pair-distribution function of water is not very sensitive 
to the value of /ijMj. 

VIII. DISCUSSION 

In the present work and in reference it has been 
shown that the inertia associated with the Kohn-Sham 
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orbitals in CPMD can affect the motion of the ions. 
There are substantial differences between the CPMD 
forces and the ground state forces at the same ionic posi- 
tions. Referencel^focussed on two materials that exhibit 
extreme behaviour. Crystalline MgO under pressure is 
an almost perfectly ionic material with a high quantum 
kinetic energy associated with valence electrons that are 
strongly bound to the oxygen ions. The /i-dependent 
error for MgO is therefore very large but can almost per- 
fectly be corrected by applying a constant mass correc- 
tion to the ions. Silicon, on the other hand, is a mate- 
rial in which valence electrons are delocalized and have 
a small quantum kinetic energy, /i-dependent errors for 
silicon are therefore extremely small. 

Here the focus has been on D2O - an important sys- 
tem in which valence electrons possess a large amount of 
quantum kinetic energy, but which is not perfectly ionic. 
It is demonstrated that /^-dependent errors are large, and 
while applying a constant mass correction to the ions im- 
proves the forces considerably, systematic errors {AF^ ) 
remain that are very large compared to the errors that 
would be present in a reasonably well-converged BOMD 
simulation. These errors are highlighted in order to un- 
derline the importance of caution when applying CPMD 
to situations in which substantial rearrangements of the 
electronic orbitals are likely to occur. The non-rigid-ion 
part of the errors, AF^, observed for crystalline D2O are 
likely to be less serious than those in liquid water, for 
example, where deviations of the electronic orbitals from 
their average structure are likely to be larger. Further- 
more, in studies of chemical reactions and phase transi- 
tions deviations from rigid-ion behavior are likely to be 
more serious and it may not be safe to assume that dy- 
namics and thermodynamics can be corrected simply by 
using mass corrections for the ions. In such situations it 
would appear safer to use a small value of ^ or to use 
BOMD. 

In the past, the quality of a CPMD simulation has 
sometimes been judged according to the degree of conser- 
vation of energy and the degree to which the adiabatic- 
decoupling condition is maintained|lOj|. In the present 
work it has been emphasized that adiabatic decouping is 
not a decoupling of the ions from the orbitals but only a 
decoupling of the ultra high-frequency part of the orbital 
motion from the lower-frequency modes of the coupled 
orbital-ion system. Perfect energy conservation and per- 
fect adiabatic decoupling have both been maintained in 
the simulations of D2O presented here, and yet the accu- 
racy of the computed forces is quite poor. 

Comparisons between BOMD and CPMD in this work 
have relied on calculating the average magnitude of sys- 
tematic errors, AF, (i.e. errors that do not average out 
on short time scales) in the forces on the ions. While 
there is no obvious alternative to using this quantity, 
and while the disparity in the magnitude of the errors 
between CPMD and BOMD appears large^H, it should 
be borne in mind that this is an extremely crude measure 
of the accuracy of a simulation. The effects of the errors 



in the forces in CPMD and BOMD are likely to be differ- 
ent. Although an analytic expression for the error in the 
forces in the case of CPMD is given by equation El the 
effect on thermodynamic properties of the part of this er- 
ror that does not reduce to a mass correction, i.e. AF'', 
is not obvious. 

In BOMD, the errors in the forces lead to a system- 
atic drift in the total energy due to a systematic bias in 
the wavcfunctions arising from their extrapolation from 
previous time steps |43i]. This energy drift can be a very 
useful judge of the accuracy of a simulation, however, a 
consequence of this is that there is also a drift in the tem- 
perature that, for long simulations, needs to be counter- 
balanced through the use of a thermostat for the ions. 
Apart from this obvious change in the temperature, the 
effects of errors in the BOMD forces on thermodynamic 
properties is unknown, in general. Recent work by Pulay 
and Fogarasi has suggested counteracting the energy drift 
by adding small corrections to the velocities of the ions at 
each time step 43] . This procedure might allow a BOMD 
simulation to be performed with a convergence criterion 
for the electronic structure that is less strict than com- 
monly used criteria as it could be both faster and with 
smaller average errors in the forces than are present in 
CPMD. A disadvantage of performing poorly converged 
BOMD simulations is that the forces on the ions would 
be discontinuous in time, and the simulations would not 
be time-reversible. In conservative CPMD simulations 
forces are always continuous and simulations are always 
time-reversible regardless of the accuracy of the forces. 

In CPMD the average kinetic energy of the ions can 
remain approximately constant throughout a simulation. 
In other words, while orbitals and ions are constantly ex- 
changing energy and momentum in CPMD, if the adi- 
abatic decoupling condition is maintained there is no 
systematic net transfer of energy between them. Ions 
both gain momentum from and lose momentum to the 
orbitals in CPMD and the physically meaningful total 
energy, which is the sum of the kinetic energy of the 
ions and the Kohn-Sham energy, fluctuates about a con- 
stant. The fluctuations in this energy are exactly mir- 
rored by the fluctuations in the fictitious kinetic energy 
because their sum is conserved. Therefore, the accuracy 
of CPMD can, in principle, be judged by the fluctuations 
in the FKE although how this should be done in practice 
is less obvious. As can be seen from figure \S\ because 
care was taken to avoid discontinuously accelerating the 
orbitals, only ionic-timescale fluctuations are visible in 
the FKE. When the rigid-ion contribution to the FKE is 
subtracted out, the FKE that remains is much smaller 
and with much smaller fluctuations, however, section IVTI 
shows that the average remaining error in the forces is 
still orders of magnitude larger than would be found in 
a reasonably well converged minimization of the orbitals 
to the electronic ground state. This is simply an illus- 
tration of the fact that small fluctuations in the energy 
(which varies only to second order in Sip) can result in 
much larger fluctuations in the force (which varies to first 
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order in Sip ). It also suggests that it may be worth ex- 
amining the effect on forces rather than energies of the 
use of thermostats to control the FKE p^[T5l |. 

Currently, ab initio molecular dynamics (AIMD) is 
hampered by two very important problems. The first 
of these is the limitations on the time scales and length 
scales that are accessible to simulations. This problem 
means that many properties and systems are currently 
beyond the reach of AIMD. For simulations that can be 
performed, this problem severely limits the precision with 
which thermodynamic properties are calculated. The sec- 
ond major problem is that the density functionals that 
are currently in use are not accurate enough for many ap- 
plications, particularly in the study of chemical reactions 
where a high level of accuracy is required on the ener- 
getics, and in systems in which electron correlation plays 
a prominent role. Both of these problems will gradually 
be reduced by technological, algorithmic and theoretical 
innovations and therefore the relative importance of the 
kinds of errors described in the present work will increase. 

It should be stressed that in the present work and in 
reference apart from clear problems with vibrational 
properties due to increased ionic masses, and changes 
in the definition of temperature, thermal pressure, and 
other quantities that depend explicitly on the mass, no 
observable thermodynamic quantity has been demon- 
strated to deviate significantly from the value that would 
be obtained in a perfectly converged BOMD simulation. 
It may be that further testing will reveal that, in practice, 
many properties are rather insensitive to /i as long as con- 
stant mass corrections are applied to the ions. Clearly, 



given the magnitude of the errors demonstrated in the 
present work, and the lack of an obvious theoretical jus- 
tification for CPMD when fi is large and the electronic 
structure deviates significantly from the rigid- ion limit, 
further testing is necessary. 



IX. CONCLUSIONS 

What has been demonstrated in the present work is 
that there are problems with some of the arguments that 
have been used in the past to theoretically justify CPMD. 
The forces in CPMD differ from the ground state forces 
even when averaged on a femtosecond time scale and even 
when constant mass corrections have been applied to the 
ions. For commonly used values of the fictitious mass, 
fi, the magnitude of these errors in the forces is large 
compared to the errors that would generally be tolerated 
in a BOMD simulation. How the errors in CPMD impact 
on thermodynamic properties is unknown. The errors 
can systematically be reduced by reducing the value of 
/i. 
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